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The effect of hydrostatic strain and of interstitial hydrogen on the elastic properties of a-iron is 
investigated using ab initio density-functional theory calculations. We find that the cubic elastic 
constants and the polycrystalline elastic moduli to a good approximation decrease linearly with 
increasing hydrogen concentration. This net strength reduction can be partitioned into a strength- 
ening electronic effect which is overcome by a softening volumetric effect. The calculated hydrogen- 
dependent elastic constants are used to determine the polycrystalline elastic moduli and anisotropic 
elastic shear moduli. For the key slip planes in a-iron, [110] and [112], we find a shear modulus 
reduction of approximately 1.6% per at.% H. 



I. INTRODUCTION 

Hydrogen degrades the performance of many alloys 
and steels by lowering the failure stress, leading to 
fracture at unpredictable loading conditionJ-32 Sev- 
eral mechanisms have been proposed to explain the H- 
cmbrittlement of iron, the main ones bein g th e hydrogen- 
enhanced decohesion (HEDE) mechanisirP^l H-vacancy 
effects^j^ and hydrogen-enhanced localised plasticity 

(HELppnni. 

In the HEDE mechanism, H weakens the cohesive 
bonds between the metal atoms, leading to failure at 
interfaces where H tends to concentrate in, suc h as 
around the tensile strain field of a crack opening} 3 -^. Va- 
cancies containing H can order thems elve s along criti- 
cal slip directions, leading to fracture^'. Within the 
HELP mechanism, the onset of plasticity with loading 
occurs at a lower stress as a result of the H-shiclding 
of repulsive interactions between dislocations^. The in- 
creased H concentration at dislocations^ effectively re- 
duces the dislocation-dislocation spacings. The resulting 
phenomenon of dislocation coalescence, and ultimately, 
crack advancement at reduced loads, in the presence of H, 
may be related to the experimentally-observed increase 
of dislocation mobility caused by EM The importance 
of increased H-concentration near dislocations and other 
low-energy trap sites such as interstitial sites was also 
indicated by a recent experimental study of intergranu- 
lar failure in steeF^. In that study, the fracture mode 
changed from ductile to brittle, as the amount of H lo- 
cated at these low-energy trap sites increased. 

One of the central challenges in avoiding hydrogen 
cmbrittlement is the interpretation of the experimen- 
tally observed effective behavioural^ that arises from 
the interplay of H solubility and diffusivity with cohe- 
sive/elastic properties, vacancies, dislocations, and other 
defects. The advantage of theroretical approaches is that 
the effects can be investigated independently, in con- 
trast to experiment. Combining highly-accurate elec- 
tronic structure calculations^^ with other methods en- 
ables access to extended time and length scales which 
are necessary for describing e.g. long-rang e stra in fields 
or small H concentrations^, kinetic effects^l2iJ am \ f or 



predicting continuum-scale propertie d 22 * 23 !. 

The goal of this study is to establish a link between 
highly-accurate ab-initio calculations and continuum 
elasticity-theory in order to explain the experimentally- 
observed influence of hydrogen on the elastic properties 
of iron (e.g. Ref.[M|) and steel (e.g. Rei.&$. To this end 
we use density-functional theory to calculate the elastic 
constants of a-Fc as a function of hydrostatic stress and 
different concentrations of interstitial hydrogen. 

In Sec. |H] we describe the details of our DFT calcula- 
tions used to calculate the elastic constants of pure Fe 
(Sec. Ill) and the modification of the elastic constants 
by H (Sec. IV). We utilize the elastic constants to cal- 
culate the effect of H on the strength parameters - bulk, 
Young's, and shear moduli - in Sec. [Viand summarise our 



findings in Sec. VI 



II. DETAILS OF THE CALCULATIONS 
A. Ab initio total energies 

The calculation of elastic properties presented here is 
based on numerical derivatives of the total energy with 
respect to strain. The large unit cells required for reach- 
ing H concentrations of a few percent makes this problem 
just within the scope of present DFT calculations. 

Our spin-polarised first-principles density functional 
theory calculations were performed using the VASP^^H 
code. W e used the projector augmented-wave 
methocP^l, with pseudopotentials considering the 3p 
electrons of Fe as valence electrons. The generalised gra- 
dient approximation (GGA) in the PW91 parametrisa- 
tiorP2 was used for the exchange-correlation functional. 
In order to verify the reliability of our conclusions, wc 
repeated some of our calculations with the Vosko-Wilk- 
Nusair^ (VWN) spin interpolation in PW91, and the 
Perdew-Burke-ErnzerhoP^ (PBE) exchange-correlation 
functional. The spin-polarised GGA has proven to 
give reliable results for the ground state 3 ^ and elastic 
properties^ for Fe. 

We used supercell geometries, a plane-wave basis with 
a cutoff of 500 eV and a T-centred k-point grid equivalent 
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to 18x18x18 for the two-atom basis bcc unit cell except 
in the case of the 128-atom unit cell where the Brillouin- 
zone sampling was equivalent to 20x20x20. We found 
these settings to be more than adequate for capturing 
the equilibrium properties, but necessary for converg- 
ing the elastic constants to within less than one percent 
error (see Sec. 
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The ions were relaxed until the 
maximum force component on each ion was less than 
0.01 eV/A while the total energies were converged to 
within 0.01 meV. The reported lattice parameter was de- 
termined from fitting the total energies to the Murnaghan 
equation of stated 

Our results for the lattice parameter (2.832 A), bulk 
modulus (194 GPa), bulk- modulus pressure derivative 
(5.42), and magnetic moment (2.17 /Ug) are in good 
agreement with other first- princi ples calculations^*^! 
and with experimental resultsP^^]. 



is applied, giving a total energy expression 

EorthoiS) = E(0) + C'VS 2 + 0[6% 

while, for C44, the strain is monoclinic 

6/2 
e= ( 6/2 

8 2 /{A~5 2 ) 



with a total energy 



1 



(5) = E{0) + -C i4 V6 z + 0[6 4 



(3) 



(4) 



(5) 



We fitted our ab initio total energies for systems sub- 
ject to nine values of strain i5 between ±3% to Eqs.|3]and 
[5] in order to obtain C and C44. 



III. ELASTIC CONSTANTS 



B. Convergence tests 



A. General description and convergence criteria 

The total energy of a solid at zero stress and equilib- 
rium volume Vq can be expanded about small strains e 



E(e) = E(0) 



yVg^Cijkieijekt + 

ijkl 



(1) 



where the indices run from 1-3. As the strain tensors are 
symmetric, the notation Cijki can be expressed in the 
two-index (Voigt) form CV, where the indices run from 
1-6. In this section we restrict ourselves to linear elastic 
behaviour, i.e. stress linear in the strain. 

For the calculation of the elastic constants in this work, 
we enforced cubic symmetry of the unit cells in our DFT 
calculations. This reflects our expectation that a ran- 
dom distribution of H would cause an effectively isotropic 
lattice distortion. In addition, the restriction to cubic 
volume elements facilitates the scale-bridging with meso- 
scopic approaches, such as e.g. finite element schemes. In 
Sec. |IVB"1 we will show that the quantitative effect of this 
cubic constraint is in fact negligible for the investigated 
unit cells. 

Within our approximation of a cubic system, the num- 
ber of independent elastic constants reduces to three: 
On, C12, and C44. Their numerical values can be deter- 
mined by applying suitable strain tensors and taking the 
second derivative of Eq. [T] with respect to the applied 
strain. The bulk modulus B = | (Cn + 2Ci 2 ) and its 
pressure derivative B'(P) are found by applying hydro- 
static strains and then carrying out a fit to the Mur- 
naghan equation of state. The values of C = C\\ — C\i 
and C44 are found by applying volume-conserving strains. 

In the case of C , an orthorhombic strain 



(2) 




To estimate the accuracy of the elastic-constant calcu- 
lations, we carried out extensive convergence tests. The 
quantities varied were the k-point density, cutoff energy, 
and strain 6 at a lattice parameter of 2.832 A . An ex- 
ample of these convergence tests is shown in Fig. [lla for 
C44, for \6\ <3% (the convergence tests for C are simi- 
lar) . The values of cutoff energy and k-point density of 
the circled point (500 eV, 18x18x18) were used in (b) 
and (c) to verify that the applied strains are within the 
linear regime. 
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FIG. 1: (a) Convergence of C44 as a function of k-point den- 
sity and energy cutoff for the two-atom bcc-Fe cell for a range 
of strain values \5\ < 3% (see Eqs. |2|5[ ). The error bars are 
deduced from the least-squares fit to Eq. [5] The circled point 
corresponds to the calculation parameters used in this study, 
(b) C44 and (c) C" obtained for different ranges of strain S. 

As a result of our extensive convergence tests, we ob- 
served less than a 2% variation in the elastic constants 
upon varying the range of 6 between 0.5-3% change in 6. 
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For the calculations presented in the remainder of this 
study, we chose a value of 3% for S, i.e. at the border 
of linear elastic behaviour, in order to avoid numerical 
instabilities and/or the need for very high precision total 
energy evaluations. 

C. Dependence of elastic constants on hydrostatic 
strain 

In order to achieve a comprehensive description of the 
influence of volume-expansion due to interstitial H, we 
also investigated the nonlinear clastic behaviour of the 
elastic constants of pure Fe as a function of hydrostatic 
strain. Equivalently, this is the effect of a volume ex- 
pansion on the linear elastic constants Cij. The modi- 
fied second-order elastic constants are given by BirctP 
in terms of the second-order and third-order elastic con- 
stants Cij and Cijk to linear order in the applied hy- 
drostatic strain. Instead of calculating the third-order 
elastic constants, we directly determined the variation in 
the second-order elastic constants as a function of small 
applied hydrostatic strain. 

Following Wallace [22] , the total energy at a volume V 
produced by strains applied at a reference volume V re f 
away from equilibrium, is modified from Eq. [ljto include 
a term first-order in strain, corresponding to hydrostatic 
stress Oij\ 

E(V,e) = E(V ref ,e = 0) + V ref J2 a ^J 

+ ^T.^ C ^i{V ref )e kl - (6) 

ijkl 

Formally, the strains are Lagrangian (second-order in dis- 
placement) strains evaluated with respect to V re * (as in 
Eq. 2.37 in [22] ) but, as we are only concerned with small 
strains, we approximate them as infinitesimal. 

Because the stress is hydrostatic, it can be written as 
crij = —PSij. Expanding Eq. [6] and allowing the strains 
represented by Eqs. [2] and |4j to be applied with respect 
to V re f , we obtain 

E ortho .(P,S) = E(P,0) + (C - P) V r ^5 2 + 0[5% (7) 



E mono .(P, 8) = E(P, 0) + - f C u - 2) yref62 + °[^ 4 ]- 

as modifications of Eqs. [3] and [5] where we have removed 
the V re f functional dependence notation from the for 
clarity. 

However, the elastic constants of Eq. [6j henceforth 
known as the energy-strain coefficients, even while mod- 
ified to be valid at V re * ', are no longer equal to the 
stress-strain coefficients. This is discussed in detail by 
Wallace 44 (Eq. 2.51) and Barron and KleirP^. Exper- 
iments mostly use either ultrasonic wave-propagation^ 



or diffraction technique s] 47 * 48 ! and obtain the pressure- 
varying elastic constants from stress-strain relations. The 

o 

stress-strain coefficients, which we call Cijki as in Rcf. 45, 
are related to the energy-strain coefficients Ciju of Eq. 
©by 

or in simplified terms 

en = Cu 

c 12 = C\2 + P 

°c i4 = C44 -\P- (10) 

We note here that the expressions [9fT0] differ from 
those foun d in works which consider Lagrangian 
strains 49 ^ but agree with others also using inifinites- 
imal strains^H^. One consequence of using the stress- 
strain coefficients instead of the energy-strain coefficients 
is that Eqs. [3] and [5] retain their form, apart from the re- 

0.00 o 

placement of the C and C44 by c =c n — c 12 and C44 

o 

respectively. We use the stress-strain coefficients Cyjy in 
order to avoid ambiguities in the definition of the Cijki 
at non-zero pressure when comparing with experiments, 
or with other computational studies. 

At each value of hydrostatic strain 77 ranging between 
±0.05, the value of 5 was varied between ±0.03 and a 
fit was made to second-order in S, as in Sec. |III A| 
These applied strain values correspond to the initial lin- 
ear regime of the high-pressure transition to the body- 
centred-tetragonal phased The elastic constants c' and 
C44 are again extracted from the coefficients of 5 2 . How- 

o o o 

ever, in order to separate Cu and c 12 from c , we need to 
take into account the change of the bulk modulus with 
pressure. The bulk modulus is modified by a term equal 
to its pressure derivative B'(P) multiplied by the pres- 
sure at the corresponding value of hydrostatic strain 77. 
The resulting variation of the stress-strain coefficients 
with hydrostatic strain follows a linear trend, as can be 
seen from Fig. [2] A linear least-squares fit to the small- 




hydrostatic strain 77 



FIG. 2: Calculated stress-strain coefficients for pure Fe as a 
function of hydrostatic strain 77 = AV/Vo with respect to the 
equilibrium volume Vo. 

strain region (|?7| < 0.02) of these curves resulted in the 
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following dependence of the stress-strain coefficients on 



hydrostatic strain: 






B(v) 


= (194- 


1075r/) GPa 


en (v) 


= (284- 


1492?7) GPa 


Cl2 (V) 


= (149- 


8877?) GPa 


C44 (77) 


= (105- 


662??) GPa 



(11) 



and the corresponding hydrostatic strain dependence 
of the energy-strain coefficients was 

6*14(77) = (284- 14927/) GPa 

C 12 (t7) = (149 - 669t/) GPa 

C u (v) = (105 - 761?/) GPa. (12) 

The difference between the stress-strain and energy- 
strain coefficients is quite small, and virtually indis- 
cernible on the scale of Fig. [2j Additionally, the use of 
another exchange-correlation functional (PW91+VWN, 
PBE) changed the absolute values of the energy-strain 
coefficients, but it had a negligible effect on the slopes of 

Pig. 21 

Previously-reported first-principles calculations of the 
stress-strain dependence of bcc-Fe (Tab.|l]) show a spread 
of approximately 10% and tend to overestimate the ex- 

o o 

perimental values for c n and c 12 while underestimating 

o 

C44. This spread may be attributed to e.g. different 
exchange-correlation functionals, applied strains, or cri- 
teria for convergence. We did an interpolation of the 
previously published data to obtain the values listed for 
the strain-dependences in the table. Experimental data 



method 


Cn(GPa) 


Ci 2 (GPa) 


c 44 (GPa) 


PAW-GGA (present) 
PAW-GGA^ 1 
LMTO-GGA^ 
PP-GGA^ 
FP-LAPVv 1 ^ 
exptP 
exptP 


284(-1492) 
271(-1228) 
303(-1282) 

289 

285 

245 

240 


149(-887) 
145(-535) 
150(-813) 

118 

139 

139 

136 


105(-662) 
101(-454) 
126(-604) 

115 

100 

122 

121 



TABLE I: Comparison of stress-strain parameters (and their 
hydrostatic strain 77 dependence in parentheses, if available) 
obtained in the present study (Eq. |ll[ | with other calculations, 
and with experimental results. 



points^HH] for the hydrostatic-strain dependencies of the 
stress-strain coefficients are displayed along with the cal- 
culated results by Sha and CoherPSI in their paper, and 
these agree well with their calculations. Our derivatives 
do not deviate appreciably compared to the other results. 

It is worth noting that the elastic moduli for a wide 
variety of phases of Fe have been found to decrease with 
applied compressive strain^. 



IV. DEPENDENCE OF ELASTIC CONSTANTS 
ON H-CONCENTRATION 

A. Simulation cells 

The elastic constants of pure Fe serve as a starting 
point for determining the influence of interstitial H atoms 
on the elastic properties. In this study, we focus on in- 
terstitial H in the tetrahedral site (Fig. [3]) because (i) 
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FIG. 3: Location of an H atom in the tetrahedral interstitial 
position on the face between two bec unit cells. Lines towards 
nearest-Fe neighbours (blue circles) are shown. 

we find this site to be 0.13 eV more stable than the oc- 
tahedral on e at z ero stress (in agreement with previous 
DFT studiei^E3) anc i (jj) we expect that the twice as 
large number of tetrahdedral sites (compared to octa- 
hedral sites) per Fe atom will dominate the mechanical 
properties at ambient temperatures. We implement the 
variation of H-concentration in the supercell approach of 
our calculations by (i) increasing the size of a supercell 
containing one H atom, or by (ii) adding a second H atom 
to the same supercell. In addition, we changed the sym- 
metry of the supercell or the position of the second H 
atom relative to the first in order to alter the ordering 
of the H atoms within the Fe host lattice. The dimen- 
sions of supercells and number of H atoms used to achieve 
various H concentrations are listed in Tab. [TTJ For each 
of these supercells we determined the elastic constants 
by assuming a cubic lattice but, by allowing for internal 
ionic relaxations, accounting in an approximate way for 
the local distortions due to H. 



B. Effect of non-cubic distortions 

The volume expansion of the Fe host lattice by H in the 
tetrahedral interstitial site introduces distortions which 
break the cubic symmetry. For the H-orientation shown 
in Fig. [3j the Fe nearest-neighbours expand radially- 
outward from H (along the red lines), resulting in a to- 
tal expansion, when projected onto the cube axes, iden- 
tical in the a and c directions, and greater than that 
along b (see e.g. Ref. I5"9"]) . This tetragonal distortion 
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at. % H 


n(H) 


k x I x m 








lxlxl 


0.8 


1 


4x4x4 


1.8 


1 


3x3x3 


2.7 


1 


3x3x2 


3.6 


2 


3x3x3 


4.0 


1 


2x2x3 


5.3 


1 


3x3x1 


5.9 


1 


2x2x2 


5.9 


1 


2x4x1 


7.7 


1 


1x2x3 


10.0 


2 


3x3x1 


11.1 


2 


2x2x2 



TABLE II: Investigated H-concentrations with corresponding 
number of H atomsi, n(H), and supercell dimensions in mul- 
tiples k, I, and m of the two-atom bcc Fe cell. For the highest 
concentration (11.1 at.%) we considered three configurations 
with H-H spacings within the supercell of 4.98 A, 3.66 A, and 
2.31 A. 



increases the number of unique elastic constants from 
three to six: C 12 = C 23 ^ Ci 3 , Cii = C33 ^ C 22 , 
and C44 = Cqq =/= C55. The variation among the no- 
longer equivalent elastic constants depends on the degree 
of tetragonal distortion, which is affected by increasing 
H-concentration and additionally may be broken depend- 
ing on the relative H positions in the supercell. 

In order to justify our use of a cubic cell, we compared 
our results with those from a tetragonal unit cell for one 
of the most distorted cases considered in our study, a 
tetragonal distortion b/a—1 of -0.4%, obtained at an H 
concentration of 11.1 atomic % and H-H spacing within 
the supercell of 4.98 A. Despite the presence of two H in 
the simulation cell, the distortion was tetragonal. The 
lattice parameters of the tetragonal unit cell were ob- 
tained from a quadratic fit over a two-dimensional grid 
of total energies. 

The elastic constants for the explicitly tetragonally- 
distorted unit cell at 11.1 at. % were calculated using the 
strain and total energy expressions given in Ref. EH1 For 
the cubic unit-cell, the elastic constants were explicitly 
calculated by straining in the different Cartesian direc- 
tions and averaging. Table |III| summarises the average 
values of the cubic and tetragonal elastic constants along 
with the associated standard errors computed from the 
spread in the values for the different orientations. Only 





cubic 


tetragonal 


Cn = 1 (Cn + C22 + C33) 

C12 = 1 (C12 + C13 + C23) 

C44 = 3 (C44 + C55 + Cee) 


240±7 
145±3 
92±2 


239±6 
147±15 
93±2 



TABLE III: Elastic parameters (GPa) for the highest- 
considered H-concentration (11.1 at.%) as obtained for a cubic 
or tetragonal (b/a-l=-0.4%) unit cell. 

C12 gave a large spread for the tetragonal cell (separately, 
the values were Ci 2 = C 23 =138 GPa, C 13 =164 GPa) 
but not for the cubic cell. The excellent agreement be- 



tween the elastic constants of the tetragonal and cubic 
cells is not surprising, as the cell volumes were found to 
be the same, which caused the variations in the elastic 
constants arising from different lattice parameters (de- 
tails in Sec. Ill C ) in the tetragonal cell to mostly av- 



erage out in the cubic cell. As a result of the excellent 
agreement with the cubic approximation for this highly- 
distorted case, we are confident that our results for all 
concentrations obtained with the cubic cell accurate. For 
the cubic unit cells with an H-concentration of greater 
or equal to 4 at.%, the elastic constants were explicitly 
calculated by straining in the different Cartesian direc- 
tions and then averaged. For lower concentrations, as 
the spread amongst the three directions was < 2%, there 
was no averaging done. Instead, the associated errors 
are from the least-squares fitting of the total energy as a 
function of strain. 



C. Single-crystal elastic constants 

The calculated variation of cubic elastic constants 
with H concentration is shown in Fig. [4j Additional 
data points at the same concentration correspond to 
differently-ordered structures (see Tab. III]). There is a 
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FIG. 4: Calculated elastic constants for bcc Fe as a func- 
tion of H concentration, along with a line of best fit to the < 
7.7% points (see text for an explanation regarding the outly- 
ing points). 

clear trend of decreasing elastic constants with increas- 
ing H concentration. We ascribe the outlying values of 
Cn and C12 for the 10 % concentration to the high stress 
associated with the particular relative orientation of the 
two H atoms, and also because of their proximity to their 
periodic images in the smallest dimension. The spread at 
11.1% is hidden in the error bars associated with the non- 
cubic distortions. At the other concentration, 5.9%, for 
which we examined different orderings, there was no sig- 
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nificant difference between the values. We constrained 
the lines of best fit shown in Fig. [4] to pass through 
the zero-concentration value, but we neglected the data 
points above 7.7% and the 0.8% value of C44. These lines 
of best fit, as a function of atomic H concentration x, are 
given by 



Cu(x) = (284 -316a;) GPa 
C 12 (x) = (149 - 72a;) GPa 
6*44(2;) = (105- 141a;) GPa. 



(13) 



Calculations with an embedded-atom potential of the 
modification of elastic moduli of Fe by H have been re- 
ported recently in Ref.[6TJ The data covered up to 6 at.% 
H and while the elastic moduli decreased with H initially, 
the effect levelled off with increasing H and generally was 
far weaker than here. 

Given the pressure-dependence of the elastic constants 
(e.g. Ref. l5"Tj) . one would expect that part of the observed 
modification of elastic properties with H concentration 
can be attributed to the volumetric effect of the inter- 
stitial H on the Fe host lattice. In our calculations, the 
volumetric effect of H is immediately apparent through 
the linear increase of supercell volume with H concentra- 
tion shown in Fig.[5j From the slope, we determined that 
each H expands the lattice by 4.5 A 3 , in excellent agree- 
ment with the value 4.4 A 3 obtained from experimental. 
This relation enables us to express the H concentration in 




atomic H concentration 



FIG. 5: Relative volume expansion as a function of H concen- 
tration where AV is the volume difference of the system with 
H compared with the volume V° of pure Fe. 

terms of a volume change corresponding to a hydrostatic 
strain rj which we can employ in our parametrisation of 
the elastic constants (Eq. 11). This correponds to a di- 
rect evaluation of the volumetric effect of H on the elastic 
constants. The volume-induced changes in elastic prop- 
erties are shown in Fig. |6j together with the total effects 
originally displayed in Fig. |4j By taking the difference 
of H-dependence and volume-dependence, we remove the 
contribution of strain of pure Fe from the elastic param- 
eters calculated at the equilibrium volume corresponding 
to each concentration of H. The resulting difference plots 
(dashed lines in Fig. [6]) show the residual effects, which 
include electronic contributions, of H at each concentra- 
tion (and implicit corresponding volume). 

The separation of solute effects from alloy elastic mod- 
uli has been also considered in Ref. [63] for H in Nb, and 
by Ref. IMl for different Fe-based binary alloys. In these 
studies, the volume effect was parametrised, by the equiv- 
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FIG. 6: Left vertical axis scale: Variation of stress-strain co- 
efficients with H (data from DFT - same as in Fig. |4| and 
volume-expected variation (parameterised from pure iron - 
see Eq. Right vertical axis scale: difference between to- 
tal and volume effect (dashed lines). 



alent of a line of best fit to the data of Fig. [5j whereas we 
determined it explicitly using the individual data points. 



V. DEPENDENCE OF STRENGTH 
PARAMETERS ON H CONCENTRATION 

A. Polycrystalline elastic moduli 

Single-crystal Fe samples with H are difficult to pre- 
pare while the measurements of the CV,- do not directly 
relate to the strength properties of the material. Most 
samples are polycrystalline, and typical measurements 
are directly related to stiffness (bulk modulus), tensile 
strength (Young's modulus), and hardness (shear modu- 
lus). In addition, microscopic simulations such as finite- 
element calculations, whose inputs consist of polycrys- 
talline averages of elastic moduli, can make direct com- 
parisons with such experiments. Therefore, we trans- 
formed our single-crystal results to poly-crystalline Fe by 

o 

using combinations of the stress-strain coefficients to 
derive various elastic moduli for describing different types 
of stress-strain responses. Shown in Fig. [7] are the bulk 
modulus, and the polycrystalline averages of the Young's 
and shear moduli, where the average was performed ac- 
cording to the expression by Hill, which is an average of 
the Voigt and Reuss bounds^. As in Fig. [6j the volu- 
metric changes in the moduli for pure Fe are also shown. 
We find a clear trend in the elastic regime towards a stiff- 
ening (B increased), tensile strengthening (E increased), 
and hardening (G increased) of the system with increas- 
ing concentration of H. Despite the overall softening of 
the material with increasing H concentration, our find- 
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ings indicate that H weakens the softening caused by the 
accompanying volume expansion. Our overall findings 
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FIG. 7: Left vertical axis scale: Variation of polycrystalline 
moduli (B- bulk, E- Young's, G-shear) with H (data from 
DFT) and volume-expected variation (parameterised from 
pure iron - see Eq. 111. Right vertical axis scale: the dif- 



ference between total and volume effect (dashed lines). 

are in reasonable agreement with the experimentally de- 
duced decrease of the shear modulus of polycrystalline 
Fe by 8% for 1 at% EP. 

The effect of H on the mechanical properties of iron has 
been controversial. A study by Matsui et a'fi' 1 found that 
in tensile tests, the flow stress is increased at low tem- 
peratures by H whereas at higher temperatures it found 
softening. The study indicated that the H-dislocation 
interaction plays an important role in determining the 
type of effect that H has on the mechanical properties. 
Specifically, whether the presence of H results in the hin- 
dering or enhancing (as in HELP) of dislocation mobility, 
seems to determine whether the effect is material harden- 
ing or softening respectively^. A very recent discussion 
of this controversy, which continues to persist, is given in 
Ref. 1681 It is important to note that the above-mentioned 
studies dealt with austenitic (fee) steel, and that similar 
experiments on ferritic (bec) steel cannot easily be con- 
ducted due to the much greater diffusivity of H in bee 
versus fee irorP^. 



B. Shear moduli in key slip planes 

Our H and volume-dependent elastic constants also en- 
able us to more closely study macroscopic failure mech- 
anisms. Therefore, we determine the H dependence of 
the shear modulus, an important quantity for describing 
the stress needed for dislocation nucleatiorfS' an d glided 
The shear modulus describes the elastic stress response to 
applied shear strain. In a single-crystal sample, the elas- 



tic moduli are anisotropic, meaning they take on different 
values when rotated to a different coordinate frame than 
the standard [100] orientation used in earlier sections. 

The shear modulus is defined by 



i^i; = 1,2,3 



(14) 



where the indices 1,2,3 denote the x, y, z axes of the co- 
ordinate system. The transformation of strains from the 
reference (unprimed) to rotated (primed) coordinate sys- 
tem can be performed using Euler angles or direction 
cosines (see e.g. Ref. 170)) . The rotated strains hence 
become e' = TeT T where T is the transformation ma- 
trix for transforming the reference coordinates x into 
the new frame x' = Tx. Simplified expressions for the 
directionally-dependent shear modulus are given in Ref. 
ITU whose notation we follow. 

Using the lines of best fit for the concentration- 
dependent elastic constants of (Fig. [4]), we are able to 
parametrise the different shear moduli as a function of 
concentration for arbitrary planes and strain directions. 
The atomic displacements during a shear distortion con- 
stitute slip, and in bec metals, the most common direc- 
tion of slip is the closest-packed (111) direction with the 
main slip planes being {110} and {112}^1. 

The two major shear moduli for cubic symmetry are 
G23 and Gyi- The G23 shear modulus describes y'z' 
shear. The diagram in the upper part of Fig. [8] shows 
the coordinate system for an arbitrary plane defining the 
x' , or 1 axis. It contains a degree of freedom, 9, defining 
the orientation of the 2 and 3 axes in the plane normal 
to the 1 axis. G12 is the shear modulus when x'y' shear 
is applied. 

We display the directionality of the shear moduli as 
polar plots in Fig. [8] for the two key slip planes of bec- 
Fe. The shear modulus G23 for pure Fe in the [HO] plane 
has a maximum for slip in the (001) and (110) direction. 
The minimum, at = 45°, does not correspond to any 
integer-multiple direction; the nearest one being (111) at 
9 « 35°, which is the expected direction of slip. For the 
[112] plane, the minimum is also the (111) direction. The 
values of G23 and G12 are the same for some directions 
but different in others. They are competing with each 
other when deciding which plane is most susceptible to 
slip and in which direction. For example, looking only at 
G23 of the [112] plane (Fig. [8]b) , it appears that both 
the (111) and (110) are equally soft, but (111) prevails 
in softness when G12 is examined. With increasing H- 
concentration, we find a nearly uniform, linear decrease 
in the shear modulus in all 9 directions. The dependence 
on 9 is weak and not visible on the scale of Fig. [8] We 
find a similar rate of decrease for all planes, amounting 
to an average over 9 of 1.6± 0.1% per atomic % H. 
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FIG. 8: Directional dependence of the shear moduli G23 and 
G12 (in GPa) for the slip planes (a) [lTO] and (b) [112] of 
bcc-Fe, for different values of H concentration. The angle 9 is 
the degree of freedom in the choice of axes in the shear plane 
(inset). 



with a tetragonally-distorted unit cell at the highest in- 
vestigated hydrogen concentration. 

From our density-functional theory calculations, we 
observe a significant linear decrease of the elastic con- 
stants Cn, C12 and C44 with increasing hydrostatic 
strain or with increasing concentration of interstitial hy- 
drogen. The volumetric part of the hydrogen dependence 
can be isolated by relating the volume dependence of the 
elastic constants to the volume expansion of the corre- 
sponding hydrogen concentration. The overall decrease 
in elastic constants is the result of two opposing contri- 
butions: the decrease in elastic constants with increasing 
volume per atom versus an increase from the electronic 
contribution. These opposing effects may help to recon- 
cile contradictory experimental findings - hardening or 
softening - under different conditions and concentrations. 

We used the elastic constants from our single-crystal 
ab-initio calculations to examine the dependence of the 
polycrystalline elastic moduli B, E, and G on the hydro- 
gen concentration and find good agreement with the few 
available measurements. The single-crystal shear mod- 
uli G\2 and G23 deduced from the ab-initio calculations 
show an isotropic decrease of approximately 1.6% per 
at.% H. This suggests that a lower yield stress (assum- 
ing the same yield strain) could be expected as a result 
of the H-lowered strength parameters in the investigated 
elastic regime. 



VI. CONCLUSIONS 

We have studied the modification of the elastic prop- 
erties of bcc Fe by hydrostatic strain and by interstitial 
hydrogen. The calculations were carried out for hydro- 
gen concentrations between 0.8 at.% and 11.1 at.% with 
simulation cells of different dimensions. Our applied con- 
straint of a cubic lattice was verified by a comparison 
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